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The R algorithm is widely used for simulating two flavours of dynamical staggered fermions. We give a simple 
proof that the algorithm converges to the desired probability distribution to within 0(5r^) errors, but show 
that the relevant expansion parameter is (5r/m)^, m being the quark mass. The Rational Hybrid Monte Carlo 
(RHMC) algorithm provides an exact (i.e., has no step size errors) alternative for simulating the square root of 
the staggered Dirac operator. We propose using it to test the validity of the R algorithm for simulations carried 
out with 5t ~ m. 



1. Introduction 

The motivation for this work is the need to de- 
cide which lattice action to use in the future, to al- 
low dynamical fermion simulations using as light 
a quark mass as possible. Recent UKQCD re- 
search has been performed using Wilson fermions, 
but the computational cost of producing ensem- 
bles is too great for mjr/mp < 0.4. One alterna- 
tive conclusion is to use two flavours of improved 
staggered fermions, since the cost appears to scale 
significantly better for small masses. 

Clearly, it is important to verify that the R 
algorithm is correct for small quark masses: this 
is the subject of this investigation. The other 
possible failings of two flavour staggered fermions 
will not be addressed here. 

2. Deriving the R Algorithm 

2.1. The $ Algorithm 

We start with the probability distribution for 
gauge field U and pseudofermion field $, 

Z 

The staggered Dirac operator is A4 with the 
pseudofermion field only defined on even sites. 

The $ algorithm is a Hybrid Molecular Dy- 
namics (HMD) algorithm for 4 flavours of stag- 
gered fermions. It iterates a composite Markov 
step, which is ergodic and has a flxed point dis- 
tribution close to the desired one. We introduce 
conjugate momenta tt in order to define a Hamil- 
tonian H , 



Z Z 

The action ScS takes the role of the potential in 
the Hamiltonian. The gauge fields U can then 
be allowed to evolve for a time r by integrating 
Hamilton's equations, using a Molecular Dynam- 
ics (MD) integration scheme. Each MD trajec- 
tory consists of a momentum refreshment heat- 
bath using Gaussian noise, a pseudofermion heat- 
bath using Gaussian noise, and finally an MD tra- 
jectory consisting of t/St steps. 

Typically a QPQ symmetric symplectic inte- 
grator is used, that is one which evolves U by 
half a step, tt by a step and finally U hy a, 
half step. This does not conserve energy, having 
SH — 0((5t^) for any trajectory length. 

For 6t > 0, the fixed point distribution of the 
MD step and the momentum refreshment heat- 
bath do not coincide. We must find the actual 
equilibrium distribution of the composite of these 
two steps. Since we discard the new momenta and 
pseudofermions after each step we consider the 
full Markov step as an update of U alone, inte- 
grating out the auxilliary fields tt and Let V{t) 
represent the evolution operator for the MD step 
V{t) : {U.tt) ^ iU",Tr") and e-(^+^^) denote 
the fixed point distribution of the full Markov 
step, where AS' measures the deviation from the 
desired distribution. This must satisfy 

g-[S(C/')+AS(C/')] 



= J dU" dn" e-^"+^^>^'' S{U' - U") 

^ g-[S(C/')+AS(C/')]^g-5(H+AS)^^^ 

with 6 : H. ^ flo [V{t) — 1] measuring the change 
in O over a trajectory (i.e., SH is the extent to 
which energy is not conserved) . We have assumed 
reversibihty V^-^ — F oV o F where F : {U, ir) i— > 
{U, — tt) and area preservation so the Jacobian is 
unity. Thus we obtain the condition 

Performing an asymptotic expansion on this con- 
dition in powers of St, knowing SH = 0(St^) for 
any trajectory length r, we deduce that 5/S.S ^ 
0{St^). We thus have AS - 0{St^) for r > St, 
hence $ is 0{St^) accurate. 

2.2. The X Algorithm 

This is very similar to <&, except the pseudo- 
fermion heatbath is performed before every sin- 
gle MD step as opposed to only before each MD 
trajectory. The proof that the leading error is 
0((5t^) follows from that of R given later. 

2.3. The Ro algorithm 

Rq begins from a completely different view- 
point. Instead of introducing pscudofcrmions to 
replace the fermionic determinant, we include the 
determinant in the action Q. The fermionic ac- 
tion is thus 

Sf = -ntrlnM^M, 

where n is the number of fermion multiplets {n = 
Nf/A for staggered fermions). When computing 
the fermionic force contribution, a noisy field es- 
timator is used for the trace. 

The ergodic composite Markov step now con- 
sists of a momentum refreshment heatbath and 
an MD trajectory consisting of = t/St steps 
with independent noise ij used for each step. 

To calculate the error in 6H we have to include 
the effect of using a noisy estimator, (S')^ — S' 
for the fermionic force. For a single noisy MD 
step we find that 

{e-'"), = - S')\il - n')Sr' + 0{Sr^). 



The coefficient of (1 — tt'^)St'^ is proportional to 
the variance of the estimated force and will only 
vanish if the force is computed exactly. Since the 
momentum average is not Gaussian after many 
leapfrog steps, the leading order term is not can- 
celled as it would be if we only did one MD 
step per trajectory. Thus over an entire trajec- 
tory SH ~ 0{St) and the leading error of Ro is 
AS ^0{St). 

2.4. The R algorithm 

The only difference between x ^^'^ Ro is that 
in the former the pseudofermion field is calcu- 
lated at the beginning of every MD step, and 
in the latter the noisy field (effectively pseudo- 
fermions) is calculated in the middle of each MD 
step. X has 0{St'^) errors for n = 1 multiplets, 
whereas Rq has errors 0{St). However for n = 
(i.e., no fermions) both are identical and have 
errors of 0{St^). We expect that the leading 
error has a linear dependence on the time the 
pscudofcrmions are refreshed and on the number 
of multiplets, so if we refresh the pseudofermion 
field at i = (1 — n)5r/2, an 0(5r^) algorithm for 
< n < 1 fermion multiplets should be obtained. 
For two flavours of staggered fermions, this means 
evaluating the pseudofermion field a quarter way 
through each MD update This leads to an 
algorithm that is neither reversible nor area pre- 
serving and so cannot be made exact (unlike the 
previous algorithms which could be made exact 
through the inclusion of an accept/reject step). 

The argument leading to Equation (|l|) may be 
generalised to give 

^^{S+S){H+AS)~tvln V-.^^ ^ X, (2) 

where S measures lack of "energy" conservation, 
S measures lack of reversibility S : i—> fl o 
[V{t)-^-FoVoF] and trln K = Indet ^^f^^ 
measures lack of area preservation. Considering 
a single step of the R algorithm, where the auxil- 
liary field x is computed at a time t = (1 — q;)(5t/2 
where a is some parameter to be determined, and 
expanding Equation (^) in St we find 

^g-(5+5)(H+AS)-trln V,^^ ^ I - ASt^ +OiST') , (3) 

where A is proportional to {n — a). If a = n the 
leading term cancels, and thus the leading error 



is 0{6t^) for the entire trajectory. Therefore, as 
claimed R is an 0((5r^) algorithm, and thus so is 
X (i.e., R with n = 1). 

3. A Source of Inaccuracy 

The staggered fermion kernel is 

where m is the fermion mass and ?7i_^ is the gauge 
field link matrix at site i in direction yu and 77^^^ 
are the staggered fermion phase factors. The St^ 
term in Equation (^) should have a coefficient 

that behaves as [{M^M)-'^-^{M''M)]^, thus 
for light modes this term could be expected to 
behave as 0{m~^). This presents no problems, 
as long as 5t is small compared to the mass. If 
6t « m, then the St expansion breaks down. For 
an exact algorithm, the accept /reject step would 
have corrected for this, however with an inex- 
act algorithm AS" ~ 0(((5r/m)^) — 0(1), so we 
would be simulating an action S+SS which differs 
from S by terms which are not small. 

When short distance observables (e.g., the pla- 
quette) are measured with St ^ m (typical of 
light fermion simulations) there is the expected 
St^ scaling, with no indication of the inaccuracy 
just highlighted. However the m'^ behaviour 
would only be expected to be true for the light- 
est fermion modes, and since bosonic observables 
do not couple strongly to these modes; we do not 
expect the behaviour to be observable here. 

Instantons correspond to zero modes of the 
Dirac operator in the massless limit, and are cru- 
cial for physical dynamical quark effects, such as 
the rj' mass. An error AS* ^ 0(1) for the lightest 
modes would thus most likely affect the instanton 
sector, and thus one may get the most interest- 
ing light dynamical quark physics wrong. Unfor- 
tunately accurate measurement of such effects is 
notoriously hard. 

4. The RHMC algorithm 

The RHMC is an exact algorithm which can 
be used to to simulate two flavours of stag- 
gered fermions. It is like HMC with two ex- 
tra ingredients: a fairly cheap but very accu- 



rate force computation and a cheap noisy ac- 
cept/reject steplDll^. 

We now write the fermion action as 

The inverse square root of the Dirac operator, can 
be approximated using an optimal Chebyshev ra- 
tional approximation. The advantage of rational 
approximations is that the error in the approxi- 
mation falls as e"/ ™ where n is the degree of 
the rational function used and m is the fermion 
mass. This accuracy is maintained over the entire 
spectrum of the Dirac operator. 

The noisy part of the accept/reject step cor- 
rects the errors in the approximation of the 
square-root, the St errors being corrected exactly. 

5. Testing the R algorithm 

We propose to test R by comparing it against 
RHMC. Initially RHMC will be used to generate 
thermalised ensembles which will then be evolved 
using R. Changes in observables and/or autocor- 
relation lengths would signal R does not get the 
correct distribution. 

Initial work has begun on the testing, although 
it is still too early to reach any conclusions. 
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